Ocean  Modelling  40  (2011)  190-198 


ELSEVIER 


Contents  lists  available  at  SciVerse  ScienceDirect 

Ocean  Modelling 

journal  homepage:  www.elsevier.com/locate/ocemod 


Nesting  a  nonhydrostatic  model  in  a  hydrostatic  model:  The  boundary  interface 

P.C.  Gallacher3’*,  D.A.  Hebert3,  M.R.  Schaferkotterb 

a  Naval  Research  Laboratory,  Ocean  Sciences  Division,  Stennis  Space,  MS  39529,  USA 
b Jacobs  Technology  Inc.,  1 020  Titan  Court,  Fort  Walton  Beach,  FL  32547,  USA 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  2  February  201 1 

Received  in  revised  form  1 1  August  201 1 

Accepted  15  August  2011 

Available  online  2  September  2011 


Keywords: 

Nonhydrostatic  model 
Boundary  conditions 
Nested  models 
Prediction 
Internal  waves 
Nonlinear  waves 


A  method  is  developed  for  adjusting  the  values  of  the  prognostic  variables  near  the  interface  between  a 
nonhydrostatic,  high  resolution  model  embedded  in  a  hydrostatic,  coarser  resolution  model.  It  incorpo¬ 
rates  a  method  of  conditioning  the  outer  domain  lateral  boundary  values  to  enforce  conservation  of  vol¬ 
ume  when  the  variables  are  interpolated  onto  the  inner  domain  grid.  This  is  accomplished  by  adjusting 
the  baroclinic  normal  velocities  at  the  open  boundaries  after  interpolation.  The  method  also  includes  a 
relaxation  scheme  which  matches  the  values  of  the  prognostic  variables  across  a  narrow  zone  near  the 
open  boundaries  of  the  submesoscale  inner  domain.  This  prevents  the  development  of  discontinuities, 
reflections  and  perimeter  currents  at  the  periphery  of  the  inner  domain.  Submesoscale  hindcasts  are  con¬ 
ducted  in  areas  of  high  Nonlinear  Internal  Wave  (NLIW)  activity.  Since  NLIWs  have  amplitudes,  0(100  m), 
which  are  not  small  relative  to  their  wavelengths,  0(100  -  1000  m)  these  hindcasts  require  a  nonhydro¬ 
static  model  in  the  inner  domain.  Since  the  hydrostatic  model  lacks  the  physics  and  resolution  to  support 
NLIWs  there  will  be  discrepancies  between  the  values  of  the  hydrostatic  and  nonhydrostatic  prognostic 
variables  at  the  boundaries.  It  is  shown  that  the  application  of  these  method  for  adjusting  the  values  of 
the  prognostic  variables  near  the  interface  allows  the  computation  of  a  nonhydrostatic,  submesoscale 
hindcast  forced  by  a  nested  hydrostatic  forecasting  system. 
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1.  Introduction 

Modeling  submesoscale  dynamics,  scales  of  O(lOm-lkm), 
requires  resolutions  that  preclude  global  modeling  and  dynamics 
that  are  unnecessary  for  global  scales.  In  order  to  model  submeso¬ 
scale  dynamics  as  part  of  a  global  forecasting  system,  as  is  the  ulti¬ 
mate  goal,  a  smaller  model  domain  with  open  boundaries  derived 
from  a  coarser,  larger  scale  model  is  required.  Thus,  some  meth¬ 
od^)  must  be  employed  to  handle  the  open  boundaries  of  the  sub¬ 
mesoscale  domain.  As  the  name  open  boundary  implies  fluid  must 
be  able  to  flow  freely  in  and  out  of  the  model  domain.  Also  scales  of 
motion  too  large  to  be  contained  in  the  submesoscale  domain  must 
correctly  influence  the  flow  in  the  domain.  The  boundaries  of  the 
submesoscale  domain  should  be  truly  transparent,  i.e.  they  should 
be  invisible  in  the  results  of  the  modeling. 

In  this  paper  we  focus  on  the  requirements  for  the  open  bound¬ 
ary  conditions  between  a  submesoscale,  nonhydrostatic  inner 
domain  and  a  mesoscale,  hydrostatic  outer  domain  for  a  one-way 
coupled  case.  The  submesoscale  resolution  allows  dynamic  fea¬ 
tures  to  be  resolved  which  are  not  captured  at  mesoscale  resolu¬ 
tion.  Furthermore,  in  order  to  model  submesoscale  dynamics 
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which  includes  Nonlinear  Internal  Waves  (NLIWs)  with  large 
amplitude  (hundreds  of  meters)  and  subinertial  eddies  and  fila¬ 
ments  a  nonhydrostatic  model  with  resolution  of  tens  to  hundreds 
of  meters  in  the  horizontal  and  meters  to  tens  of  meters  in  the  ver¬ 
tical  is  required.  The  difference  in  the  physics  between  the  nonhy¬ 
drostatic  model  and  the  hydrostatic  model  will  also  lead  to 
differences  in  the  solutions  that  exacerbate  the  problems  of  inter¬ 
facing  the  nested  models.  When  a  mismatch  of  the  solutions  occurs 
at  the  interface  between  the  domains  spurious  reflections  and 
erroneous  perimeter  currents  along  the  boundaries  (Mason  et  al., 
2010)  can  result.  Also  incorrect  geostrophic  currents  can  be  set 
up  along  the  open  boundaries  by  erroneous  density  gradients 
across  the  boundaries. 

The  open  boundary  conditions  (OBCs)  must  be  effective  at  pro¬ 
viding  information  to  the  inner  domain  from  the  outer  domain 
and  at  propagating  information  out  of  the  inner  domain.  Both 
requirements  must  be  accomplished  transparently.  In  one-way  cou¬ 
pling,  which  is  the  focus  of  this  paper,  the  information  propagated 
out  of  the  inner  domain  does  not  affect  the  outer  domain.  Because 
of  the  differences  in  physics  and  resolution  there  will  be  inherent 
differences  between  the  mesoscale  solution  and  the  submesoscale 
solution.  Prognostic  variables  from  the  coarse  grid  must  be  interpo¬ 
lated  onto  the  fine  grid  in  a  manner  that  conserves  basic  properties 
and  scales.  This  interpolation  can  be  problematic  since  the  most 
common  interpolation  schemes  are  not  conservative,  e.g.  bilinear 
and  nearest  neighbor  (Guo  et  al.,  2003  and  Accadia  et  al.,  2003). 
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Formulating  and  applying  OBCs  to  limited  area  models  is  chal¬ 
lenging  and  complex.  This  has  been  an  active  area  of  research  for 
some  time,  see  reviews  by  Palma  and  Matano  (1998)  and  Davies 
(1983).  OBC  methods  can  be  divided  into  two  categories:  (1)  adap¬ 
tive  and  (2)  consistent.  Adaptive  methods  differentiate  between 
incoming  and  outgoing  fluxes  and  apply  different  OBCs  values 
based  on  the  direction  of  the  flux.  Consistent  OBCs  apply  the  same 
value  regardless  of  the  direction  of  the  boundary  flux. 

The  advantage  of  using  adaptive  OBCs  is  that  different  OBCs  can 
be  applied  for  each  case.  This  allows  the  response  to  be  tailored  to 
the  characteristics  and  scales  of  the  outer  domain  or  inner  domain 
solution  as  appropriate  (Marchesiello,  2001).  Disadvantages  of 
adaptive  OBCs  include:  determining  which  points  qualify  as  inflow 
and  which  qualify  as  outflow,  handling  transitions  between  adja¬ 
cent  inflow  and  outflow  areas  and  handling  strong  tangential 
flows.  Adaptive  OBCs  are  further  complicated  by  the  fact  that  any 
open  boundary  point  could  simultaneously  be  an  inflow  point  for 
some  scales  and  an  outflow  point  for  others. 

On  the  other  hand,  consistent  OBCs  apply  the  same  method  to 
all  open  boundary  points.  These  methods  are  typically  a  form  of 
wave  radiation,  absorption  or  relaxation  scheme.  Problems  with 
radiation  schemes  stem  from  the  fact  that  multiple  waves  with  dif¬ 
ferent  phase  speeds  will  exist  in  realistic  cases;  whereas  the  radi¬ 
ation  schemes  use  only  one  phase  speed.  Absorption  schemes 
typically  have  a  region  where  the  viscosity  is  artificially  increased 
to  damp  the  solution  near  the  boundaries.  Relaxation  schemes  in¬ 
volve  a  region  near  the  boundaries  in  which  the  inner  and  outer 
values  are  smoothly  merged.  One  issue  with  relaxation  and  absorb¬ 
ing  layers  is  that  the  equations  governing  the  motion  in  the  bound¬ 
ary  layer  are  not  the  correct  equations  for  any  approximation  of  the 
physics  so  spurious  solutions  can  result.  See  Blayo  and  Debreu 
(2005)  for  a  discussion  of  OBCs  based  on  characteristics.  In  this  pa¬ 
per  we  focus  on  methods  for  consistent  OBCs  and  the  particulars  of 
the  methods  are  discussed  below. 

For  simulations  involving  idealized  or  climatological  case  stud¬ 
ies  of  processes  the  large  scale  flow  can  be  prescribed  as  constant 
or  a  simple  function  of  time  and  the  OBCs  merge  these  values  with 
the  flow  field  from  the  submesoscale  domain.  For  predictions 
whether  in  real-time,  i.e.  forecasts,  or  offline,  i.e.  hindcasts,  the 
problem  is  more  complicated.  In  these  cases  the  OBCs  must  be  able 
to  merge  the  multiple  time  and  space  scales  of  the  larger  scale 
dynamics  and  the  multiple  time  and  space  scales  of  the  resolved 
submesoscale  dynamics. 

The  remainder  of  this  paper  is  divided  as  follows:  Section  2 
which  discusses  the  boundary  conditioning  methods  we  have  em¬ 
ployed,  Section  3  which  describes  the  model  and  the  domain  for 
the  experiments,  Section  4  which  describes  the  results  of  several 
experiments  demonstrating  these  boundary  conditioning  methods, 
and  Section  5  which  summarizes  the  results  and  gives  some 
conclusions. 


2.  Boundary  conditioning 


near  the  interface.  The  FRS  allows  outer  domain  values  to  propa¬ 
gate  into  the  inner  domain  with  minimal  distortion. 


2.3.  Transport  Correction  Scheme 


The  outer  domain  values  are  interpolated  to  the  inner  domain 
grid  at  points  along  the  lateral  boundaries  using  bilinear  interpola¬ 
tion.  i;  is  also  interpolated  from  the  coarse  grid  to  the  fine  grid. 
Additionally  since  the  outer  domain  is  a  part  of  an  operational  sys¬ 
tem  the  values  may  have  been  interpolated  onto  different  grids  for 
archiving  or  analysis  purposes  prior  to  the  interpolation  onto  the 
inner  domain  grid.  Each  interpolation  can  generate  errors  resulting 
in  a  loss  of  volume  conservation.  For  both  the  transport  and  t /,  and 
particularly  the  former,  the  interpolation  is  not  conservative  (Guo 
et  al.,  2003).  Furthermore,  the  fact  that  the  transport  and  i(  are 
interpolated  separately  can  lead  to  an  imbalance  between  the 
fields  (Shulman  et  al.,  2002).  Small  errors  in  the  normal  velocities 
through  the  sides  can  cause  significant  changes  in  transport  into 
the  inner  domain  and  this  is  reflected  in  large  deviations  in  )/  over 
time. 

Volume  conservation  requires  that  the  net  transport  through 
the  lateral  boundaries  is  equal  to  the  temporal  change  of  rj  in  the 
inner  domain.  However,  the  geometry  of  ocean  domains,  i.e.  small 
aspect  ratio  (h/L)  guarantees  that  the  change  in  17  will  be  small. 
Therefore,  we  correct  the  transport  through  the  domain  by  setting 
the  average  transport  to  zero  for  every  time  interval  between  the 
outer  domain  boundary  values  (one  hour  in  this  case). 

The  average  transport  through  the  lateral  boundaries  can  be 
specified  as: 

St,  =  Zi  *  T]  J  jr  Jo  u„{b,z)dzh(b)dbdt  (1) 


where  b  are  the  points  along  the  open  boundary,  r,  un  is  the  normal 
baroclinic  velocity  at  b,  h  is  the  depth  at  b  and  t,  and  t2  are  succes¬ 
sive  outer  domain  data  time  points.  The  change  in  )/  averaged  over 
the  inner  domain  is  given  by 


dJfjdS 

dt  JsdS 


(2) 


where  S  is  the  surface  area  of  the  inner  domain.  If  volume  is  con¬ 
served  then  fT]  =  Tjt.  However,  as  we  have  discussed  above  interpo¬ 
lation  and  other  factors  prevent  this  equality  from  holding.  Thus  we 
have 


St,  =  n t  +  k  ~  et  (3) 

where  is  the  error  in  the  transport  and  it  is  much  larger  than  TJt. 
This  imbalance  is  corrected  by  adjusting  un  uniformly  at  every  grid 
point  along  the  open  boundaries  at  every  inner  domain  timestep 
based  e£.  Thus 

u*(b,z)  =  un(h,z)-p|p.h  (4) 


We  have  extended  and  improved  two  OBC  methods  which  to¬ 
gether  correct  the  problems  of  conservation  and  matching  at  the 
interface  between  a  hydrostatic  outer  domain  and  a  nonhydrostat¬ 
ic  inner  domain.  The  first  method  is  the  Transport  Correction 
Scheme  (TCS).  OBCs  are  formulated  to  conserve  volume  by  match¬ 
ing  the  change  in  sea  surface  elevation,  )/,  and  the  total  transport 
through  the  lateral  boundaries  and  distributing  the  difference  into 
the  baroclinic  transports.  The  second  method  is  the  Flow  Relaxa¬ 
tion  Scheme  (FRS).  In  the  Flow  Relaxation  Scheme  OBCs  are  con¬ 
structed  to  prevent  reflections  and  perimeter  currents  which  can 
be  generated  when  flows  do  not  exit  the  interior  domain  correctly. 
This  is  related  to  a  mismatch  of  coarse  and  fine  resolution  values 


is  the  corrected  baroclinic  normal  velocity  from  the  outer  domain 
interpolated  onto  the  open  boundary  of  the  inner  domain  where  h 
is  the  inward  directed  normal  vector  at  the  boundary. 

2.2.  Flow  Relaxation  Scheme 

Differences  in  resolution,  physics,  and  numerics  will  cause 
hydrodynamic  features  to  be  different  in  scale  and  dynamics  be¬ 
tween  the  outer  and  inner  domains.  No  errors  that  arise  due  to  vari¬ 
ations  between  these  values  can  be  allowed  to  propagate  into  the 
inner  domain.  To  accomplish  this  a  relaxation  zone  is  included  in 
the  inner  domain  adjacent  to  the  lateral  boundaries.  The  relaxation 
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zone  is  an  area  where  the  fluxes  of  conserved  baroclinic  quantities 
from  the  coarse  and  the  fine  grids  are  matched. 

Relaxing  the  solutions  can  be  viewed  as  two  separate  problems. 
First  allowing  outwardly  directed  fluxes  to  exit  smoothly  and  with¬ 
out  reflection  or  the  generation  of  artificial  perimeter  currents.  Sec¬ 
ond  allowing  the  fluxes  from  the  outer  domain  to  enter  and  affect 
the  inner  solution  without  creating  a  discontinuity  or  being  exces¬ 
sively  damped.  Each  problem  is  addressed  below. 

The  first  problem  has  been  addressed  by  radiation  or  diffusion 
of  conserved  quantities.  Radiation  OBCs  such  as  Flather  et  al. 
(1976),  or  Orlanski  (1976)  assume  that  a  wave,  usually  barotropic, 
propagates  the  conserved  quantities  out  of  the  inner  domain.  The 
problem  with  Flather  and  Orlanski  OBCs  is  that  there  is  only  one 
phase  speed  so  multiple  waves  and  non-wave  phenomena,  such 
as  eddies  and  filaments  are  not  propagated  out  of  the  domain 
effectively.  Nevertheless  the  Flather  OBC  has  had  some  success  in 
realistic  three  dimensional  hydrostatic  models.  In  fact  the  Flather 
boundary  condition  insures  near  conservation  for  mass  and  energy 
through  the  open  boundary  (Blayo  and  Debreu,  2005).  For  example, 
it  is  used  when  nesting  Navy  Coastal  Ocean  Model  (NCOM) 
domains  and  is  sufficient  to  maintain  volume  conservation  in  the 
inner  domain  (C.  Rowley,  personal  communication)  without 
explicitly  enforcing  conservation  of  baroclinic  values  at  the  lateral 
boundaries.  However,  this  did  not  work  for  the  nonhydrostatic 
submesoscale  domain  we  are  implementing.  Thus  we  had  to  con¬ 
struct  a  new  set  of  OBCs. 

Diffusive  techniques  employ  a  form  of  sponge  or  absorbing 
layer  where  the  viscosity  and  diffusivity  are  artificially  increased, 
damping  the  inner  domain  conserved  quantities  to  zero  as  the 
boundary  is  approached.  The  damping  needs  to  be  smooth  and 
complete  such  that  there  is  no  reflection  back  into  the  domain 
and  no  generation  of  spurious  modes  in  the  relaxation  zone. 

The  second  problem  requires  a  technique  for  matching  the  solu¬ 
tions  across  the  boundary.  A  method  to  accomplish  this  is  the  Flow 
Relaxation  Scheme  (FRS).  Several  schemes  were  reviewed  by  Da¬ 
vies  (1983)  and  were  further  described  by  Martinsen  and  Engedahl 
(1987).  The  FRS  method  was  expanded  by  Lavelle  and  Thacker 
(2008),  hereafter  referred  to  as  LT08.  Various  geometric  forms  for 
the  relaxation  coefficient  have  been  proposed  by  Jensen  (1998) 
and  Modave  et  al.  (2010)  (hereafter  referred  to  as  MDD10).  All  of 
these  formulations  were  for  barotropic  OBCs,  here  we  apply  the 
method  to  the  baroclinic  OBCs. 

The  basic  formulation  for  the  FRS  is 


D X 
Dt 


RHS-o(n)X 


(5) 


where  /  is  any  of  the  baroclinic  prognostic  variable,  e.g.  tempera¬ 
ture,  T,  salinity,  S,  eastward  velocity,  u,  or  northward  velocity,  v. 
RHS  represents  all  the  terms  on  the  right  hand  side  of  the  Eulerian 
equation  describing  the  total  rate  of  change,  D/Df,  of  the  prognostic 
variable  x  and  cr(n)  is  the  relaxation  coefficient  in  the  relaxation 
zone.  Note  that  the  relaxation  zone  is  part  of  the  inner  domain  solu¬ 
tion  with  the  addition  of  the  relaxation  term.  In  the  relaxation  zone 
the  coordinate  normal  to  the  boundary  is  n.  The  width  of  the  relax¬ 
ation  zone  is  &  =  NAn  with  n  =  0  at  the  interior  edge  of  the  relaxation 
zone  (Fig.  1).  Note  that  a  can  also  be  a  function  of  the  tangential 
coordinate  in  the  relaxation  zone  as  in  the  “simple”  relaxation 
scheme  of  LT08.  However,  as  noted  by  LT08,  variation  of  a  tangen¬ 
tially  is  only  necessary  in  the  case  of  poorly  defined  outer  domain 
values.  In  our  case  those  values  are  well  defined  since  we  are  using 
a  nested  model  solution  rather  than  climatology. 

MDD10  defines  a  characteristic  length  scale  for  the  damping  of 
features,  i.e.  the  width  of  the  relaxation  zone,  in  terms  of  a  baro¬ 
tropic  gravity  wave  speed  and  the  damping  coefficient  as 


*  10'3 


Fig.  1.  Schematic  of  the  relaxation  zone  at  western  boundary  showing  value  of  o (n) 
where  n  is  the  coordinate  normal  to  the  boundary  and  width  of  relaxation  zone 
(3  =  NAn,  N=  10)  gridpoints).  Similar  relaxation  zones  exist  inside  each  of  the  four 
lateral  boundaries  between  the  inner  domain  and  the  outer  domain. 


(6) 


We  generalize  (6)  to 


(7) 


where  un  is  the  nominal  normal  velocity  of  the  dominant  baroclinic 
feature  moving  through  the  relaxation  zone.  If  l„  is  the  nominal 
length  scale  of  the  features  to  be  modified  in  the  relaxation  zone, 
then  the  width  of  the  relaxation  zone  should  be  approximately 
<5  =  NAn  >  ln  where  N  is  order  10.  Then  we  suggest  that  the  optimal 
absorption  coefficient  scales  like 


MDD10  and  several  other  authors  argue  that  the  relaxation  zone  is 
more  effective  if  a  varies  spatially  across  the  relaxation  zone. 
MDD10  compares  several  geometric  forms  for  a  including  a  polyno¬ 
mial  form 


ff(n)  =  (JmQ  (9) 

where  a  >  0  and  am  is  the  value  at  the  outer  edge  of  the  relaxation 
zone.  When  a.  =  2  this  reduces  to  the  quadratic  form  of  LT08  for 
their  “simple”  sponge  layer.  We  tested  that  form  and  a  hyperbolic 
tangent  form, 

o{n)  =  (Jm  tanh  Q ,  (10) 

which  was  suggested  by  Martinsen  and  Engedahl  (1987)  and  Jensen 
(1998).  We  found  that  the  hyperbolic  tangent  worked  best. 

Based  on  the  above  discussion  and  some  experimentation  the 
relaxation  zone  in  our  model  is  10  grid  points  wide.  The  optimal 
value  of  am  is  0.004  for  our  case  based  on  tests  described  below. 
The  spatial  variation  of  a  is  given  by  (10)  (see  Fig.  1,  which  also 
shows  the  grid  convention  in  the  relaxation  zone).  There  is  a  relax¬ 
ation  zone  adjacent  to  each  of  the  four  boundaries. 


3.  Models  and  domains 

The  nested  configuration  (Fig.  2)  used  for  these  experiments 
consists  of  an  inner  domain  using  the  MIT  model  in  a  2°E  x  1°N 
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Fig.  2.  Luzon  Straits  1/64°  (LZS64)  RELO  NCOM  domain  with  the  Dongsha  Plateau  1/200“  (DSP200)  MIT  domain  shown  by  the  red  rectangle  (Dongsha  Island  is  the  red  dot  in 
the  southwestern  corner  of  the  DSP200  domain). 


domain  centered  at  (21°N,  117.75°E)  with  1/200°  (approximately 
500  m)  horizontal  resolution  which  we  name  the  DSP200  (Dongsha 
Plateau  1/200°)  domain.  This  is  nested  inside  an  outer  domain 
using  a  8°E  x  6°N  RELO  NCOM  (Coelho  et  al.,  2009)  domain  at  1/ 
64°  (approximately  2  km)  horizontal  resolution  which  is  named 
the  LZS64  (Luzon  Straits  1/64°)  domain  (Chao  et  al„  2007).  The 
LZS64  domain  obtains  its  initial  conditions  and  boundary  condi¬ 
tions  from  archived  data.  The  data  was  archived  during  forecasts 
by  a  59°E  x  45°N  regional  NCOM  domain  at  1/16°  (approximately 
6-9.5  km)  resolution  designated  the  EAS16NFS  (East  Asian  Seas 
1/16°  Nowcast/Forecast  System)  (Riedlinger  et  al.,  2006).  The 
EAS16NFS  is  part  of  the  1/8°  (nominally  14  km  resolution)  Global 
NCOM  Nowcast/Forecast  system  (Barron  et  al.,  2006). 


Table  1 

Model  parameters. 


Latitude  (x) 

Longitude  (y) 

z 

Nominal  grid  spacing 

1/200° 

1/200° 

0.75- 

(500  m) 

(500  m) 

270  m 

Number  of  grid  points 

348 

200 

35 

Viscosity  (mV1) 

0.24 

0.24 

7.0  x  10~4 

Temperature  diffusivity 
(mV1) 

3.4 

3.4 

1.0  x  10~4 

Salinity  diffusivity  (mV1) 
Time  step  (s) 

3.4  x  10~2 

5 

3.4x1 0-2 

1.0  x  10~6 

The  MITgcm  nonhydrostatic  ocean  model  (Marshall  et  al., 
1997a,  b)  is  an  open  source  code  which  is  available  to  the  commu¬ 
nity,  (http://mitgcm.org/).  It  is  a  finite  volume  model  with  orthog¬ 
onal  curvilinear  coordinates  in  the  horizontal  and  variable  cell 
thickness  in  the  vertical.  Several  options  are  available  for  advection 
and  diffusion  calculations.  In  nonhydrostatic  mode  a  three  dimen¬ 
sional  equation  must  be  solved  for  the  pressure.  This  is  done  using 
a  preconditioned  conjugate  gradient  solver.  The  model  includes  a 
complete  treatment  of  the  Coriolis  force  and  a  partial  cell  treat¬ 
ment  of  topography  (Adcroft  et  al.,  1997).  The  model  uses  a  flexible 
domain  decomposition  method  for  efficient  parallel  processing, 
allowing  it  to  be  run  on  modern  supercomputers.  The  model 
parameters  for  this  study  are  given  in  Table  1. 

The  global  NCOM  nowcast/forecast  system  is  forced  with  the 
Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS) 
(Hogan  and  Rosmond,  1991)  and  with  a  regional  forecast  system, 
the  Coupled  Ocean  Atmosphere  Prediction  System  (COAMPS)  (Ho- 
dur,  1997).  Also  equilibrium  tides  are  obtained  from  the  OSU  (Eg¬ 
bert  et  al.,  1994)  global  tidal  database.  Ten  constituents  are 
included:  I<1,  01,  PI,  Ql,  M2,  S2,  N2,  K2,  MF,  MM.  This  includes 
both  height  and  velocity  values  of  the  tides.  The  bathymetry  used 
for  the  system  is  the  2  min  NRL  Digital  Bathymetric  DataBase  (NRL 
DBDB2,  http://www7320.nrlssc.navy.mil/DBDB2_WWW/),  a  global 
2-min  ocean  bathymetry  data  base  which  is  based  on  the  NAVO 
global  5-min  DBDBV  bathymetry  and  is  enhanced  with  other 
higher  resolution  bathymetries  and  coastlines.  The  model  also 
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assimilates  satellite  altimeter  data  and  Multi-Channel  Sea  Surface 
Temperature  (MCSST)  from  satellite-derived  AVHRR  to  further  im¬ 
prove  accuracy. 

4.  Results 

Submesoscale,  nonhydrostatic  hindcasts  were  conducted  in  the 
DSP200  domain  for  30  days  starting  on  April  1,  2005  at  0000Z.  The 
hindcasts  were  initialized  using  archived  data  from  LZS64  and  the 
boundary  conditions  were  also  from  LZS64  data  archived  at  one 
hour  intervals.  The  initial  data  included  baroclinic  tides,  long 
wavelength  internal  waves  and  fronts  that  develop  into  NLIWs  in 
the  nonhydrostatic  domain. 

4.1.  Transport  Correction  Scheme 

Initially  we  coupled  the  domains  using  a  Flather  boundary  con¬ 
dition  on  the  barotropic  flow  and  no  explicit  conservation  scheme 
for  the  baroclinic  flow  field.  This  is  a  standard  practice  in  coupled 
ocean  models  and  conserves  mass  sufficiently  for  most  cases 
(Blayo  and  Debreu,  2005)  and  (C.  Rowley,  personal  communica¬ 
tion).  However,  in  our  case  this  produced  significant  phase  errors 
in  f/  which  caused  us  to  examine  the  sea  surface  elevation  varia¬ 
tions  in  more  detail.  The  change  in  spatial  mean  sea  surface  eleva¬ 
tion  in  the  inner  domain,  which  results  from  the  transport 
based  on  (1)  is 

1?  =  Qr,  (11) 

The  values  of  i/f,  which  decrease  to  a  minimum  of  30  m  and  then 
increase  to  a  maximum  of  nearly  18  m,  are  unrealistically  large 
and  indicate  erroneously  large  changes  in  volume.  This  is  typical 
of  boundary  conditions  with  non-conservative  transport  as  dis¬ 
cussed  in  Section  2  above. 

The  TCS  described  in  Section  2.1  was  applied  to  the  normal 
velocities  interpolated  from  the  outer  domain  at  the  lateral  bound¬ 
aries  of  the  inner  domain.  The  resulting  JT  rj?dt  is  shown  in  Fig.  3. 
The  values  are  much  smaller  than  the  uncorrected  values  indicat¬ 
ing  that  the  corrected  transport  achieves  a  much  better  approxi¬ 
mation  of  volume  conservation.  Large  oscillations  or  large  trends 
in  i/f  can  be  viewed  as  an  indicator  of  the  error  in  the  volume 
transport  which  resulted  from  the  interpolation  of  the  outer  do¬ 
main  variables  onto  the  inner  domain  grid.  As  itf  becomes  smaller 


it  exhibits  more  realistic  values  and  temporal  variability  such  as 
the  tidal  variations,  which  can  be  seen  in  Fig.  3. 

Without  the  TCS  the  values  of  t/P  were  -30  m  to  20  m.  With  the 
TCS  we  obtain  values  for  p?  of  0.6  to  1.4  m.  Point  measurements  of 
)/  in  the  northern  South  China  Sea  yield  values  of  0.3  m  to  0.7  m 
(Beardsley  et  al.,  2004)  and  (Lien  et  al„  2005)  and  tidal  models  give 
0.6-1. 0  m  (Jan  et  al.,  2007).  The  values  of  i may  still  be  slightly 
large;  however,  the  they  are  in  much  better  agreement  with  point 
observations  than  before.  The  TCS  also  reduces  erroneous  trends  in 
thermocline  displacement  and  volume  averaged  temperature  in 
the  model  results. 


4.2.  Flow  Relaxation  Scheme 

Ideally  the  open  boundaries  of  a  nested  model  should  be  trans¬ 
parent  to  both  incoming  and  outgoing  flows.  However,  some  sub¬ 
mesoscale  features,  e.g.  eddies  and  NLIWs  are  sometimes  reflected 
at  the  open  boundaries  where  they  should  propagate  out  of  the  in¬ 
ner  domain  (Fig.  4,  top  row).  This  is  due  to  the  mismatch  of  prop¬ 
erties,  such  as  propagation  speed  and  amplitude,  between  the 
inner  domain  representation  of  the  features  and  the  outer  domain 
representation  of  the  features.  The  increased  resolution  in  the  in¬ 
ner  domain  allows  features  to  develop  with  smaller  scales,  larger 
amplitudes  and  different  phase  speeds  than  the  outer  domain  rep¬ 
resentation  of  the  same  features.  Also  features  will  develop  that  are 
not  resolved  in  the  outer  domain.  This  is  due  to  improved  dynam¬ 
ical  resolution  and  lower  numerical  viscosity.  Furthermore,  the 
nonhydrostatic  physics  includes  a  more  complete  representation 
of  the  vertical  momentum  resulting  in  more  robust  vertical  motion 
in  the  inner  domain  and  additional  features. 

In  addition  to  reflections,  mismatches  between  the  inner  do¬ 
main  and  outer  domain  solutions  can  cause  the  formation  of  cur¬ 
rents  flowing  along  the  perimeter  of  the  inner  domain.  We  find 
some  evidence  of  perimeter  currents  in  our  numerical  experiments 
usually  associated  with  reflections.  However,  the  reflections  are 
more  abundant  and  much  more  disruptive.  Also  we  find  that  elim¬ 
inating  the  reflections  also  eliminates  the  perimeter  currents  so  we 
focus  on  eliminating  reflections  from  the  boundaries. 

The  FRS  discussed  in  Section  2  was  employed  to  reduce  or  elim¬ 
inate  these  problems.  The  FRS  matches  the  outer  domain  solution 
and  the  inner  domain  solution  in  a  relaxation  zone  adjacent  to  the 
interface  between  the  domains  (Fig.  1).  Experiments  designed  to 
determine  the  optimal  FRS  were  compared  graphically  and  using 


Fig.  3.  Average  sea  surface  height,  t],  after  application  of  the  Transport  Correction  Scheme  (TCS).  The  smaller  range  of  i;  indicates  improved  volume  conservation.  The  tidal 
signal,  which  is  a  strong  component  of  the  forcing  in  this  area,  is  now  clearly  visible. 
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Fig.  4.  Temperature  at  150  m.  Top  row  is  the  experiment  with  no  FRS.  Bottom  row  is  the  FRS  experiment.  In  panels  A  and  A'  a  NLIW  is  clearly  visible  in  the  middle  of  the 
panel.  In  B  and  B'  the  waves  propagate  toward  the  western  boundary.  In  C  and  C  the  waves  arrive  at  the  western  boundary.  In  D  the  wave  has  reflected  into  the  domain.  In  D' 
there  is  no  reflected  wave. 


the  total  kinetic  energy  averaged  over  the  inner  domain  for  the 
entire  experiment, 

m  =  (u2  +  v2  +  w2)dVdt  (12) 

as  a  quantitative  metric.  KE  is  chosen  as  a  metric  because  waves 
artificially  reflected  from  the  boundaries  will  increase  the  KE  in 
the  inner  domain.  This  is  explored  further  below. 

The  reflection  of  NLIWs  at  the  lateral  boundaries  of  the  inner 
domain  can  be  seen  most  clearly  at  the  western  boundary.  The 
hydrostatic  representation  of  the  NLIWs  generally  enter  the  inner 
domain  through  the  eastern  boundary.  They  are  transformed  by 
the  nonhydrostatic  physics  and  improved  resolution  in  the  inner 
domain  and  by  interactions  with  the  topography.  More  NLIWs 
are  locally  generated  and  these  interact  with  the  remotely  gener¬ 
ated  NLIWs  as  they  propagate  across  the  inner  domain  in  a  mainly 
westerly  direction.  Thus  the  mismatches  between  the  low  resolu¬ 
tion,  hydrostatic  solution  of  the  outer  domain  and  the  high  resolu¬ 
tion,  nonhydrostatic  solution  in  the  inner  domain  are  particularly 
large  at  the  western  boundary.  As  the  waves  attempt  to  exit  the  in¬ 
ner  domain  this  mismatch  of  solutions  creates  reflections  which 
are  very  obvious  in  animations  of  the  inner  domain  solutions  of 
the  experiments  without  the  FRS. 

One  of  the  reflection  events  is  shown  using  the  temperature 
field  at  150  m  on  the  western  side  of  the  domain  (Fig.  4).  At  hour 
67  (Fig.  4A  and  A')  a  NLIW  front  is  clearly  seen  in  the  middle  of 
the  panel  in  both  the  FRS  experiment  (Fig.  4A')  and  the  experiment 
without  FRS  (Fig.  4A).  Without  the  FRS  the  NLIW  front  approaches 


the  western  boundary  and  is  reflected  (Fig.  4A-D).  The  front  is 
most  obvious  north  of  Dongsha  Island.  At  hour  69  (Fig.  4B  and 
B')  the  front  is  approaching  the  western  boundary.  At  hour  71 
(Fig.  4C)  it  begins  to  intersect  the  western  boundary  and  at  hours 
73  it  can  be  seen  reflecting  from  the  western  boundary  and  prop¬ 
agating  back  into  the  inner  domain  in  the  case  without  the  FRS 
(Fig.  4D).  The  wave  then  travels  back  into  the  central  region  of 
the  inner  domain  and  significantly  impacts  the  solution. 

In  the  experiment  with  the  FRS  (Fig.  4A'  through  D')  no  part  of  the 
front  is  reflected  from  the  western  boundary.  At  hour  69  (Fig.  4B'),  a 
wave  front,  similar  to  that  seen  in  Fig.  4B,  approaches  the  western 
boundary.  At  hour  71  (Fig.  AC)  it  starts  to  intersect  the  western 
boundary.  No  reflected  front  is  visible  at  hour  73  (Fig.  4D').  The  solu¬ 
tions  with  the  FRS  and  without  the  FRS  are  not  identical.  This  is  ex¬ 
pected  since  there  is  less  reflected  energy  in  the  FRS  experiment  than 
in  the  experiment  without  the  FRS.  Also,  since  the  energy  levels  are 
different,  the  details  of  the  solution  in  the  inner  domain  will  not  be 
the  same  and  the  propagation  speeds  of  the  NLIWs  and  Submeso- 
scale  Features  (SMFs)  will  also  be  different. 

In  addition  to  the  qualitative  visual  assessment  of  the  FRS,  we 
have  determined  that  the  magnitude  of  KE  can  be  used  as  a  quan¬ 
titative  measure.  If  NLIWs  and  SMFs  are  artificially  reflected  at  the 
open  boundaries  then  the  magnitude  of  KE  should  be  larger  in  the 
inner  domain.  Values  of  KE  are  listed  in  Table  2,  where  the  largest 
I<E  was  found  for  the  experiment  with  no  FRS.  This  confirms  our 
expectations  that  spurious  reflections  at  the  boundaries  will  lead 
to  higher  energy  levels. 

Experiments  were  then  conducted  to  minimize  I<E  using  the 
quadratic  and  hyperbolic  tangent  profiles,  (9)  and  (10),  and  for 
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Table  2 

Average  kinetic  energy  in  inner  domain  as  a  function  of  maximum  relaxation 
coefficient  and  relaxation  profile. 


Relaxation  profile 

Maximum  relaxation  coefficient, 

(7m  x  102 

K£(x  100) 

No  relaxation  region 

N/A 

9.32 

Quadratic 

4 

8.89 

Hyperbolic  tangent 

4 

8.80 

Hyperbolic  tangent 

2 

8.79 

Hyperbolic  tangent 

1 

8.78 

Hyperbolic  tangent 

0.4 

8.77 

Hyperbolic  tangent 

0.2 

8.78 

Hyperbolic  tangent 

0.1 

8.82 

Hyperbolic  tangent 

0.02 

8.93 

several  values  of  the  maximum  relaxation  coefficient,  erm.  The  re¬ 
sults  are  shown  in  Table  2jmd  Fig.  5.  From  Table  2  it  can  be  seen 
that  in  all  FRS  cases  the  I<E  was  lower  than  in  the  case  with  no 
FRS.  The  best  results  for  the  FRS  are  obtained  with  a  small  but  finite 
am  and  a  hyperbolic  tangent  profile. 

There  does  exist  a  am  which  yields  a  minimum  I<E  as  can  be  seen 
in  Fig.  5  and  Table  2.  At  one  extreme  am  cannot  become  arbitrarily 
small  since  as  <rm  approaches  zero  the  case  with  no  FRS  is  recov¬ 
ered,  see  (10).  However,  neither  can  am  become  arbitrarily  large. 
This  is  a  result  of  the  discrete  nature  of  the  model.  As  am  becomes 
larger  the  spatial  variation  of  the  amplitude  of  the  damped  wave  in 
the  relaxation  zone  becomes  too  large  to  resolve  with  a  given  An. 
This  is  different  from  the  continuous  case  where  am  can  become 
arbitrarily  large  since  dn  can  be  infinitesimal  (MDD10). 

The  spatial  variation  of  a  also  matters.  If  a  changes  too  slowly 
across  the  relaxation  zone  the  wave  will  not  be  completely  damped 
before  it  passes  through  the  relaxation  zone  and  will  be  reflected 
off  the  outer  domain  edge  of  the  relaxation  zone.  If  ct  changes 
too  quickly  the  wave  will  be  reflected  near  the  inner  domain  edge 
of  the  relaxation  zone.  This  is  equivalent  to  moving  the  boundary 
close  to  the  inner  edge  of  the  relaxation  zone.  A  small  correction 
which  changes  gradually  at  the  interior  edge  of  the  relaxation  zone 
and  more  rapidly  as  the  outer  domain  is  approached  provides  the 
best  result.  The  hyperbolic  tangent  profile  fits  this  criteria  better 
than  the  quadratic  profile.  This  is  in  agreement  with  the  results 
of  MDD10  and  shows  the  importance  of  a  spatially  varying  relaxa¬ 
tion  coefficient. 

Although  it  is  important  to  reduce  or  eliminate  the  reflection  of 
NLlWs  and  other  SMFs  at  the  inner  domain  boundaries,  it  is  also 
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Fig.  6.  (A)  Temperature  from  inner  domain  solution  and  (B)  outer  domain  solution 
cropped  to  inner  domain  area  (bottom  panel)  at  95  m  for  hour  221. 

important  to  transmit  information  from  the  outer  domain  into 
the  inner  domain  without  altering  the  signals  in  either  magnitude 
or  phase.  To  check  this  we  compare  the  inner  domain  solution  to 
the  outer  domain  solution  cropped  to  the  inner  domain  area 
(Fig.  6).  The  details  of  the  two  snapshots  are  different;  however, 
the  large  scale  patterns  are  very  similar.  This  indicates  that  the 
information  is  being  passed  into  the  inner  domain  has  maintained 
the  low  wavenumber  characteristics  of  the  outer  domain  solution. 

To  quantify  this  we  compared  wavenumber  temperature  spec¬ 
tra  for  the  outer  domain  solution  cropped  to  the  inner  domain  area 
to  that  for  the  inner  domain  solution.  First  we  took  data  along 
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Fig.  5.  Average  horizontal  kinetic  energy  (KE)  as  a  function  of  the  maximum  absorption  parameter  (<rm). 
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Fig.  7.  Temperature  spectra  from  the  inner  domain  solution  (dotted  line)  and  the  outer  domain  solution  (solid  line)  cropped  to  the  inner  domain.  Spectra  were  calculated 
along  longitudes  then  averaged  across  latitudes. 


longitudinal  grid  lines,  using  grid  lines  starting  10  grid  points  into 
the  inner  domain  from  the  southern  boundary  and  ending  10  grid 
points  into  the  inner  domain  from  the  northern  boundary.  We  cal¬ 
culated  the  spectrum  along  each  of  these  lines  then  averaged  the 
spectra  over  all  lines.  The  resulting  spectra  are  very  similar  up  to 
the  wavenumber  cutoff  for  the  outer  domain  (Fig.  7).  The  spectra 
from  the  inner  domain  continues  to  higher  wavenumbers.  This 
demonstrates  that  the  low  wavenumber  energy  has  been  passed 
from  the  outer  domain  solution  to  the  inner  domain  solution  with 
no  appreciable  alteration. 


5.  Conclusions 

A  nonhydrostatic  model  has  been  successfully  embedded  in  a 
hydrostatic  model  as  part  of  a  prototype  forecast/hindcast  system. 
The  key  is  a  carefully  crafted  treatment  of  the  prognostic  variables 
near  the  open  boundaries.  Both  a  TCS,  to  conserve  volume,  and  a 
FRS  to  merge  the  outer  domain  and  inner  domain  solutions  near 
the  boundary  are  required  for  successful  nesting.  Hindcasts  in 
the  inner  domain  show  high  wavenumber,  large  amplitude  fea¬ 
tures  which  are  not  present  in  the  outer  domain  due  to  missing 
nonhydrostatic  physics  and  resolution.  The  features  in  the  inner 
domain  agree  with  point  observations.  Obtaining  area  averaged 
time  series  data  to  verify  is  a  topic  for  future  work.  Verification 
of  additional  model  data,  e.g.  T,  S,  u,  fand  w,  with  in  situ  and  remo¬ 
tely  sensed  data  are  in  progress  and  will  be  reported  in  another 
paper. 

Domains  must  be  nested  Matryoshka  (Russian  Doll)  style  to 
achieve  the  resolution  needed  to  model  submesoscale  dynamics 
as  part  of  a  global  system.  To  accomplish  realistic  hindcasts  the 
submesoscale  domain  must  be  nested  in  a  mesoscale  domain, 
which  is  nested  in  a  regional  domain,  which  is  nested  inside  a  glo¬ 
bal  domain.  Through  the  combination  of  real-time  forecast  systems 
and  off-line  hindcast  systems  we  can  achieve  a  global  to  submeso¬ 
scale  modeling  system. 


Open  boundaries  naturally  occur  when  modeling  submesoscale 
dynamics  since  the  domains  of  interest  are  limited  by  the  small 
grid  spacing  needed  to  resolve  the  dynamics.  The  boundaries  must 
provide  information  about  the  larger  scale  outer  domain  dynamics 
on  the  inner  domain  grid  without  introducing  erroneous  trends. 
Also  the  mismatches  of  the  solutions  in  the  inner  and  outer  do¬ 
mains  must  be  adjusted  at  the  boundaries  to  prevent  contamina¬ 
tion  of  the  inner  domain  solution.  These  mismatches,  caused  by 
differences  in  resolution,  physics  and  numerics,  manifest  at  the 
boundaries  of  the  inner  domain  as  reflections  of  momentum  and 
energy  and  as  spurious  perimeter  currents. 

OBCs  from  the  outer  domain  must  be  mapped  onto  the  inner 
domain  grid  and  must  be  adjusted  to  conserve  volume.  We  used 
a  TCS  to  adjust  the  baroclinic  normal  velocities  such  that  the 
change  in  ij  is  balanced  by  the  total  barotropic  transport  through 
the  open  boundaries.  The  baroclinic  normal  velocities  are  weighted 
by  the  lateral  surface  area  of  the  grid  cell  adjacent  to  the  open 
boundaries  to  conserve  total  transport.  This  method  improves 
model  performance.  Techniques  for  improving  this  method  and 
gathering  observations  for  comparison  are  topics  of  ongoing  and 
future  work. 

Dynamic  features  generated  or  transformed  in  the  inner  domain 
can  be  quite  different  from  those  in  the  outer  domain.  In  fact  some 
dynamic  features  in  the  inner  domain  will  not  be  present  in  the 
outer  domain.  This  can  be  due  to  increased  resolution  in  the  inner 
domain,  in  which  case  the  features  are  not  resolved  in  the  outer 
domain.  It  can  be  due  to  differences  in  hydrostatic  versus  nonhy¬ 
drostatic  physics,  in  which  case  the  feature  would  not  exist  in 
the  outer  domain  even  if  the  resolution  in  the  two  domains  was 
the  same.  Finally  this  can  be  due  to  differences  in  the  numerics 
or  the  numerical  parameters  between  the  two  domains,  in  which 
case  the  features  might  be  damped  out  or  dissipated  in  the  outer 
domain.  In  any  case  the  differences  will  be  manifested  as  discrep¬ 
ancies  in  the  prognostic  variables  near  the  open  boundaries.  This 
can  lead  to  reflections,  perimeter  currents  and  other  spurious  fea¬ 
tures  of  the  solution  on  the  inner  domain.  These  spurious  features 
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cannot  be  allowed  to  propagate  into  the  inner  domain  and  contam¬ 
inate  the  results.  A  solution  is  to  construct  a  relaxation  zone  along 
the  perimeter  of  the  inner  domain  where  these  spurious  features 
are  removed  from  the  inner  domain  solution  using  a  FRS. 

It  is  desirable  to  have  the  relaxation  zone  as  narrow  as  possible 
to  maximize  the  area  available  for  the  inner  domain  solution.  At 
the  same  time  the  spurious  features  must  be  damped  well  enough 
that  they  do  not  significantly  impact  the  solution  in  the  inner  do¬ 
main  or  the  boundary  values  from  the  outer  domain.  The  spatial 
characteristics  and  maximum  value  of  relaxation  coefficient  are 
both  important  to  the  effectiveness  of  FRS  in  the  relaxation  zone. 
We  have  shown  that  a  relatively  small  relaxation  coefficient  with 
a  hyperbolic  tangent  spatial  distribution  that  varies  slowly  near 
the  inner  domain  side  of  the  relaxation  zone  and  varies  rapidly 
near  the  outer  domain  side  of  the  relaxation  zone  yields  the  best 
results.  This  is  in  agreement  with  previous  findings  by  other 
authors. 

Nonhydrostatic  physics  is  required  to  obtain  the  correct  dynam¬ 
ics  for  the  submesoscale  domain.  This  increases  the  differences  be¬ 
tween  the  inner  domain  and  outer  domain  solutions  which 
complicates  the  task  of  reconciling  the  two  solutions  near  the 
boundaries.  This  necessitated  a  careful  treatment  of  both  the  TCS 
and  the  FRS. 
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